W problemie tym poszukujemy skalarnego pola ciśnienia
\( {\cal R}^3 \ni (x,y,z)\times [0,T] \rightarrow u(x,y,z;t) \in {\cal R } \)
w obszarze niejednorodnym. Problem ten związany jest z zagadnieniem wydobycia ropy naftowej. Zakładamy iż obszar w którym modelujemy nasz przepływ składa się materiału niejednorodnego reprezentującego złoża ropy rozmieszczone w szczelinach warstw geologicznych. Zakładamy iż wykonano odwiery w których umieszczono pompę pompującą ciecz pod wysokim ciśnieniem, której celem jest wypompowanie złoża ropy naftowej w szczególności w kierunku pomp ssących umieszczonych w oddzielnych odwiertach.
Równanie modelujące zmiany ciśnienia na wskutek działania pomp
\( \frac{\partial u}{\partial t } - \nabla \cdot (\kappa \exp(\mu u)\,\nabla u) = f \)
co po szczegółowym rozpisaniu operatorów różniczkowych
\( \frac{\partial u(x,y,z;t)}{\partial t } - \frac{\partial }{\partial x} \left(\kappa(x,y,z) \exp(\mu(x,y,z) u(x,y,z;t)) \frac{\partial u(x,y,z;t)}{\partial x} \right) + \frac{\partial }{\partial y} \left(\kappa(x,y,z; u(x,y,z;t)) \exp(\mu u(x,y,z;t)) \frac{\partial u(x,y,z;t)}{\partial y } \right) = f(x,y,z;t) \)
W równaniu powyższym funkcja \( \kappa(x,y,z) \) określa przepuszczalność obszaru. Funkcja ta przedstawiona jest na rysunku Rys. 1. Zakres wartości zmienia się od 0 do 1000. Wartość 0 oznacza brak przepuszczalności, wartość 1000 oznacza bardzo wysoką przepuszczalność (np. szczelinę w skale w której znajduje się ropa naftowa).

Współczynnik \( \mu=1000 \) jest parametrem modelu.
Zauważmy że równanie przepływu nieliniowego jest bardzo podobne do równania transportu ciepła. Jedyna różnica to nieliniowy współczynnik
\( \kappa(x,y,z) \exp(\mu(x,y,z) u(x,y,z;t) \)
który oznacza że prędkość zmian ciśnienia w danym obszarze zmienia się proporcjonalnie do współczynnika przepuszczalności oraz do aktualnej wartości ciśnienia, i dzieje się to w sposób nieliniowy, określony funkcją eksponencjalną. Współczynnik ten jest przyczyną nieliniowości naszego równania. Gdyby zostawić ten człon na lewej stronie równania, wówczas równanie to wymagało by zabiegów linearyzacji na przykład za pomocą algorytmu Newtona-Raphsona. W przypadku problemów nieliniowych, możliwe jest jednak zastosowanie metody explicite która daje możliwość przesunięcia członu nieliniowego na prawą stronę równania.
Pochodną po czasie przybliżamy za pomocą metody różnic skończonych.
\( u_{t+1}-u_t - dt* \nabla \cdot (\kappa \exp(\mu u)\,\nabla u) = dt*f \)
Funkcja prawej strony modeluje pompy pompujące i pompy ssące.
\( f(x, t) = \textcolor{blue}{\sum_{p \in P} \phi\left(\left\|x_p - x\right\|\right)} - \textcolor{red}{\sum_{s \in S} u(x, t) \phi\left(\left\|x_s - x\right\|\right) } \)
gdzie \( \textcolor{blue}{\sum_{p \in P} \phi\left(\left\|x_p - x\right\|\right) } \) to pompy pompujące (dostarczające niezerowego ciśnienia) oraz \( -\textcolor{red}{\sum_{s \in S} \phi\left(\left\|x_p - x\right\|\right) } \) to pompy ssące ("dostarczające" ujemnego ciśnienia). Suma przebiega po wszystkich położeniach pomp pompujących P i ssących S w naszym modelowanym układzie.
Ciśnienie wymuszane przez pompa spada z odległością od środka pompy zgodnie ze wzorem, zilustrowanym również na rysunku Rys. 2

\( \phi(t) = \begin{cases} \left(\frac{t}{r} - 1\right)^2 \left(\frac{t}{r} + 1\right)^2 & \text{for } t \leq r \\ 0 & \text{for } t > r \\ \end{cases} \).

W celu przeprowadzenia symulacji za pomocą izogeometrycznej metody elementów skończonych wprowadzamy dyskretyzację za pomocą funkcji bazowych B-spline
\( \{B^x_{i,2}(x)B^y_{j,2}(y)B^z_{k,2}(z)\}_{i=1,...,N_x;j=1,...,N_y;k=1,...,N_z } \)
Będą one stosowane do aproksymacja symulowanego pola skalarnego ciśnienia w aktualnej chwili czasowej
\( u(x,y,z;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y;k=1,...,N_z} u^{t+1}_{i,j,k } B^x_i(x)B^y_j(y)B^z_k(z) \)
Podobnie do testowania użyjemy funkcji bazowych B-spline:
\( \{ B^x_l(x)B^y_m(y)B^z_n(z) \}_{l=1,...,N_x;m=1,...,N_y;n=1,...,N_z } \)
Nasze równanie przemnażamy przez funkcję testującą i odcałkowujemy po obszarze trójwymiarowej symulacji (całkowanie oznaczamy nawiasem reprezentującym iloczyn skalarny w przestrzeni L2)
\( (u_{t+1},w)=(u_t,w) - dt ( \nabla \cdot (\kappa \exp(\mu u)\,\nabla u),w)+dt (f,w) \)
Możemy odcałkować przez części człon na prawej stronie, zakładając że nasz obszar jest tak duży iż na brzegu zmiany ciśnienia wynoszą zero (wtedy zniknie nam całka brzegowa wynikająca z całkowania przez części)
\( (u_{t+1},w)=(u_t,w) + dt ( \kappa \exp(\mu u)\,\nabla u,\nabla w)+dt (f,w) \).

Podobnie jak w ogólnym przypadku metody explicite, na lewej stronie mamy teraz macierz masową którą można zfaktoryzować w czasie liniowym, a na prawej stronie mamy naszą "bitmapę" która opisuje nam stan w chwili poprzedniej plus zmiany w przeciągu kroku czasowego dt wynikające z fizyki przepływu nieliniowego w ośrodku niejednorodnym oraz z modyfikacji pól ciśnienia w okolicy pomp pompujących i ssących.